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1 0 BACKGROUND OF THE INVENTION 

X-ray computed tomography (CT) has been very successful in imaging internal 
structures of the human body. It has provided an accurate, micro-resolution and 
real-time medical imaging tool for clinical use. However, the use of X-rays has several 
disadvantages. The contrast is low for certain kinds of tumors, such as early stage 

15 breast cancer. The misdiagnose rate for X-ray mammography is high and X-rays are 
mutagenic. In recent years, imaging techniques based on optical photons have attracted 
significant interest. In the spectral region between 700 and 900 nm, called the 
therapeutic window, photons do not give rise to mutagenic effects, and they can 
penetrate deeply into tissues, due to the weak absorption of light at these wavelengths. 
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Sensitivity to optical contrast is high and spectroscopic information can be obtained. 
Further, contrast can be enhanced by injecting exogenous dyes which target tumor cells. 
Thus, the information provided by optical photons can complement that of X-ray CT, 
and perhaps provide an alternative diagnostic tool for detecting tumors and other 
5 abnormalities inside the body. 

SUMMARY OF THE INVENTION 

The major obstacle to optical imaging is the high turbidity of human tissues. 
Photons injected into tissue undergo multiple scattering events before they are detected. 
To extract spatial information, one has to disentangle the effects of scattering. Various 

1 0 optical imaging approaches have been developed to deal with this problem. Each of 
these approaches has advantages and disadvantages. Imaging techniques utilizing 
ballistic photons, which are very good for imaging up to one millimeter inside the 
tissue, do not work in the case of deep tissue imaging. Several image reconstruction 
schemes have been devised to deconvolve turbidity and improve spatial resolution. 

1 5 Finite element methods and finite difference methods have been developed in the time 
domain and the frequency domain, respectively. Filtered backprojection has been 
applied to CW imaging. Weighting factors for the contribution of individual voxels to 
the measurement, first proposed in X-ray CT, have also been calculated using Monte 
Carlo simulations and employed in the inverse model. The present invention relates to 

20 the use of series expansion methods to the optical regime. Using an early time detection 
approach, three dimensional images of a tissue can be reconstructed by taking into 
account the effects of turbidity. 

The problem can be understood by comparing the propagation of optical photons 
and X-rays in human tissue. As it is well known, X-rays traversing the body (Figure la) 

25 travel along straight lines. In contrast, due to scattering, the propagation of optical 

photons has a three dimensional spread. The distribution of optical photon paths can be 
visualized as a tube connecting the source and the detector (Figure lb). The width of 
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the cross section of the tube varies according to the time at which arriving photons are 
collected. 

Several features of the photon path distribution can be noted: 
First, in the time-of-flight limit this tube reduces to a straight line, and the 
5 problem is reduced to that of conventional X-ray CT. However the signal produced by 
these so-called ballistic photons is extremely weak, and essentially non-detectable for 
thick tissues; Second, for early detection times (a few hundred ps after the time of 
flight), the tube is still very narrow and the photon paths are well defined. The 
transmitted signals are significant, and spatial information is still highly preserved; and 
10 Third, at late detection times (several ns), the tube can completely fill the organ of 

interest, and spatial resolution is significantly reduced. This regime is referred to as the 
diffusive limit. 

By replacing the straight line paths of X-rays with the narrow tubes of the early 
time photon path distribution, an optical CT procedure can be employed that is a 

1 5 modification of that used in X-ray CT. 

In the optical regime, the early arriving photons are analogous to X-ray photons. 
They undergo a smaller number of scattering events in comparison with highly diffusive 
photons, and thus preserve a significant amount of spatial information. Yet signal levels 
for early arriving photons can be relatively high. Measurements taken using early time 

20 detection have higher resolution compared to those obtained with continuous wave 
(CW) and frequency-domain techniques. Sharp images can be reconstructed using the 
concept of photon path density. As is well known, the diffusion approximation solution 
does not well describe the early arriving photons. The predictions of the diffusion 
approximation with a point spread function up to the t-t 0 -600 ps time window are 

25 inadequate for correct image reconstruction. A point spread function (PSF) that takes 
into account causality produces much better results. The correct form of the point 
spread function for early time plays an essential role for image reconstruction. 

The image reconstruction method of the present invention is based on the use of 
a series expansion method in the optical regime, where scattering is dominant and the 
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distribution of photon paths between source and detector must be taken into account. 
To accomplish this, a PSF is used to generate a weighting function matrix. Previously, 
weighting functions have been discussed and calculated using the diffusion 
approximation, the microscopic Beer-Lambert law, and Monte Carlo simulations. 
5 However, the use of the PSF provides guidance into the choice of weighting functions. 
The physical interpretation is clearer in terms of the PSF, and different theories about 
photon migration can be tested because they predict different PSF's. Furthermore, as 
shown in Eq. (12) below, there are actually two sets of weighting functions, one for 
scattering contrast and one for absorption contrast. For example, tumors in breast tissue 

10 exhibit both absorption and scattering contrast. The early portion of the photon 

migration curve is more sensitive to the scattering contrast than the absorption contrast. 

The resolution of optical CT is restricted by several factors, such as the effects of 
scattering and the underdetermined nature of the reconstruction procedures. 
Additionally, the total number of projections and measurements can be increased, and a 

1 5 fan-beam geometry can be used to improve data collection efficiency. Fiberoptic 
systems can be used for delivery and/or collection of light from the tissue of a patient 
under examination. Refined time-domain photon migration instruments, implemented 
with a computer using reconstruction programs can provide optical CT images with 
high quality in the breast, the brain and elsewhere in the body. 

20 BRIEF DESCRIPTION OF THE DRAWINGS 

The foregoing and other objects, features and advantages of the invention will be 
apparent from the following more particular description of preferred embodiments of 
the invention, as illustrated in the accompanying drawings in which like reference 
characters refer to the same parts throughout the different views. The drawings are not 

25 necessarily to scale, emphasis instead being placed upon illustrating the principles of the 
invention. 

Figures 1 A and IB are schematic diagrams employing an algebraic 
reconstruction technique for Optical CT in which the sample under study is divided into 
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N x N x N vowels and the absorption distribution is represented by the average 
absorption within each voxel 

Figures 2 A and 2B illustrate each voxel being assigned a weighting factor 
including for the X-ray, in which the voxel on the trajectory has 100% contribution, 
5 while off the trajectory has 0% contribution and for the optical photons, the weighting 
factor is determined by the photon path distribution, respectively. 

Figures 3 A and 3B are schematic diagrams of systems for performing optical 
computed tomography in accordance with the invention. 

Figure 4 illustrates an example of the dimension of the scattering medium and 
1 0 the scanning geometry. 

Figures 5 A and 5D are contour maps of the intensities for 1 -object (a)-(b) and 
2-objects (c)-(d) with time windows of: (a)-(b) At = 534 ps, width 50ps; (c)-(d) At = 607 
ps, width 50 ps. 

Figures 6A-6D are point spread functions in which the source detector are at (0.- 
15 3.2.0) cm, and (0.3.2.0) cm, respectively, and four combinations are given, (A) the side 
view of the z=0 plane, calculated with the diffusion approximation, (B) the side view of 
the Z=0 plane, calculated with the causality correction: (C) the top view of the y=0 
plane, calculated with the diffusion approximation; (D) the top view of the y=0 plane, 
calculated with the causality correction. 
20 Figures 7A-7D are reconstructed images with one embedded object including 

(A) reconstruction with direct X-ray algorithm; (B) reconstruction with diffusion 
approximation; (C) reconstruction with causality correction; (D) the exact configuration. 

Figures 8A-8D are images of two embedded objects including (A) 
Reconstruction with direct X-ray algorithm; (B) reconstruction with diffusion 
25 approximation; (C) reconstruction with causality correction; (D) the exact configuration. 
Reconstruction with the direct X-ray algorithm does not resolve the two embedded 
objects. On the other hand, the reconstructions with diffusion approximation overdo the 
de-convolution and result in smaller images. The smaller object is invisible from the 
3D view. 
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Figures 9A-9D show the contour map diagram of four slices of the reconstructed 
image in Figure 8C where the slices are (a) Z = 10mm; (b) Z = -6mm; (c) Z = 2mm; and 
(d) Z = 10mm. 

Figure 10 illustrates a process sequence of a preferred embodiment of the 
5 invention. 

DETAILED DESCRIPTION OF THE INVENTION 

The representation of the attenuation of the X-ray intensity through the human 
body is well established. Assume the human tissue under study has an attenuation 
distribution ji* (r 1 ) for a monochromatic X-ray. Then the transmitted intensity of the 

1 0 X-ray beam across the body is : 

a) 



I(r) = I Q exp 



0 a 



where the integral is along the straight line connecting the source at r 0 and the detector 
at r. Equation (1) can be rewritten as: 
15 (2) 



r 

-[ZnZ(r)-In/ 0 ] = J^(r'). 



The above line integral is often referred to as the Radon transform. 

Unlike X-rays, near infrared light in human tissue undergoes strong scattering 
and does not follow straight-line paths. Optical photons undergo multiple scattering 
20 events before they are detected or absorbed by the tissue. The distribution of photon 
paths in a uniform scattering medium has been studied using various approaches. It has 
been established that the distribution of the photon paths is narrow at early detection 
times, but spreads out as time increases. 
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Generally, the presence of tumors creates optical heterogeneity which appears as 
a local variation of the scattering and the absorption properties of the tissue. For early 
stage tumors, these changes can be considered as small perturbations. 

The propagation of near infrared photons in human tissue is well described by 
the transport equation. For a medium with scattering distribution \i s (r), and absorption 

distribution \x a (r), the radiance Z(r,s,t ) satisfies: 



(3) 



-^^+5-VI(r,5,0 = -[^(r) + ^ s (r)]l(r,5,0 

c ct L 



+fi, (r)JU(r,s, t)f r (5-5' )d£l'+Q(r , s, t) 



4 it 



and the normalized differential scattering cross-section / (s * s' )> often referred to as 
1 0 the phase function, satisfies : 

1 I f (S-s')dQ'=l 
4% 



(4) 



In Eq. (3), c is the speed of light in the medium, and Q(r, s,t) is the source term. 
For a perturbative system, the distribution of absorption, scattering and phase function 
1 5 have small variations of the form: 

(5) 

Va( T ) = f t a +S Ma( r )> 

fi s (r)/ r (s • j ' ) = juj(s • j* ) + Sju s (r,s ■ s' ), 



with SM s (rJ-S') = Sju s (r)f(s-s') + ju s [f r (S-S')-f(s-S')]. 
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In Eq. 5, \i a and ju s are the average absorption and scattering coefficients, respectively; 

8[i a (r)(« // a ) and 8 j\^f s (r)(« juj are the perturbations caused by the tumors. 

Generally, the local variation of the phase function f T (s • §') from the global phase 
function/(s * §') may not be small for some angles, but Eq. (4) requires that such 
5 variations cancel each other after integrating over all solid angles. Thus treat 

5ju s (rj- s) as a small perturbation. 

Under variations given by Eq. (5), Eq. (3) can be solved using perturbation 
theory. For the case of a point source, 

(6) 

10 fi(r,5,0 = ^-^(r-r 0 )^-0» 

Eq. (3) is the equation for the Green's function G (r,t| r 0 1 0 ), which has an 
expansion 

(7) 

G i (r,*|r 0 ,< 0 ) = Gi 0) (r,/ 0 |r o ,/ o ) + Gi 1) (r^°|r o ,r o )+... 
15 with Gj- 0 \r,t\r 0 ,t 0 ) being the normal Green's function in the absence of perturbation 
and Gf\r,t\ r 0 ,t 0 ) the first order correction due to the perturbation. Substituting Eqs. 
(5)-(7) into Eq. (3) and keeping terms only up to the first order, we have: 

(8) 

c at 

-fi,llGW(r,t\r 0 ,t a )f(.i-i')da'= ^-S{r - r 0 )<5(t-t 0 ) 

20 and 
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(9) 



1 dGj» (r ;iV 0 ) +i . VG (i )(r>f | v ) + ( + ^)GO) (r ,f|r 0 f 0 ) 
-^JjGj, 1) (r,r|r 0> g/(i.l')dn' 

= -[<$u(a)(r) + Sju s (r)]Gf } (r, f|r 0 ,t 0 ) + ll G£ 0) (r, f|r 0 , t 0 )tfi s (r,s- s)dQ' 



Equation (9) can be solved through the adjoint equation of Eq. (8), which is: 



1 dG™(r 9 t\r 0 t 0 ) 

c dt 0 



■^V 0 Gf(r^|r o g + (^ +// s )Gf (r^|r 0 / 0 ) 



(10) 



-ju s J J GP (r,t\r 0 ,t 0 )f(s-s')da=-t-S(r - r 0 )5(t - 1 0 ). 



It can be shown that the solution to Eq. (9) satisfies: 

(11) 

^IirfQGjV^,0 = -f'^ fdV[ < y// 1 (r0 + ^,(r')]JJdGi o) (r,f|r o ,f o )Gi 0) (r' s r|r o ,f o ) 

J t 0 J 4^r 4^ 



Under the assumption of uniqueness of the solution to the transport equation, 
10 Eq. (1 1) can be reduced to: 
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(12) 

-]-G?\r,t\r 0 ,t 0 ) = -f df \^d^)Gf\r,t\r\f)Gf\r\f\r () ,Q 

47T to 

-f df Jd 3 r'^ s (r')Gf (r,t\r\f)G^(r\f\r 0 ,t 0 ) 

to 

+ f df f </Vj ldQ'Gl°\rj\r\nG^\r\t%J 0 )bJi s (r\s ■?) 

"To J 4x 

-f df§dA'-sGP {rj\r\f)Gf\r\f\r,J o y 



The first term on the right hand side of Eq. (12) is the correction due to 
absorption variations, the second and the third terms contain the correction due to 
5 scattering variations; the third term also includes an extra correction due to phase 

function variations. The last term in Eq. (12) is the surface integral and is related to the 
boundary condition. For boundary conditions commonly used, such as the 
zero-boundary condition in our system described below, the surface integral contributes 
at higher order and can be discarded. In order to simplify the formulation, define the 
1 0 point spread function as: 

(13) 

PSF(r,t;r';r 0 ,t 0 ) = 4nf dfGl°\r,t\r\f)Gi Q \r\f\r 0 ,t 0 ). 



Therefore, Eq. (7) and Eq. (12) can be rewritten as: 



15 



,t n )-Gi°\r,t 



0' 0 



Vo } 



= fd 3 r^ a (r')-PSF(r,t^;r 0 t 0 ). 



(14) 



Note the remarkable similarity between Eq. (2) and Eq. (14). Physically, 
PSF (r J ; r ' ; r 0 ,t 0 ) is the probability that a photon is inj ect ed at the source point r 0 at 
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Note the remarkable similarity between Eq. (2) and Eq. (14). Physically, 
PSF(r 9 t; r 1 ; r 0 ,t 0 ) is the probability that a photon is injected at the source point r 0 at 

time t 0> passes through the field point r' before time t, and is collected by the detector at 
point r at time t. It represents the photon path distribution at time t between the light 
5 source and the detector. In the ballistic limit, the photon path distribution in the 
transverse direction (i.e., direction perpendicular to r-r 0 ) shrinks to a S-function: 

(15) 

PSF(r,f;r';r 0 ,0) Hr " r ° |/c > 5 x (r-r Q ). 

Thus, the volume integral on the r.h.s, of Eq. (14) reduces to a line integral along 
10 r-r 0 , and Eq. (2) and Eq. (14) are essentially identical. 

It is important to note that the derivation of Eq. (14) does not employ the 

assumptions of the diffusion approximation. The Green's function G| 0) (r ; /|rV) in 

Eq. (13) can be calculated using various models. Three examples are the path integral 
solutions, the random walk solution, and the conventional diffusion approximation 

1 5 solution. Further details regarding time gated imaging methods can be found in U.S. 
patent No. 5,919,140, the entire contents of which is incorporated herein by reference. 
In fact, these reconstructions show that the conventional diffusion solution is inadequate 
for imaging with early arriving photons. It predicts a point spread function which is too 
wide and renders images which are too small. A solution satisfying causality is required. 

20 For the purpose of comparison, consider the diffusion approximation: 

(16) 

Gi 0) (r^|r 0 ,/ 0 )s^G (0) (r^|r 0 ^ 0 )--^-— 5-VG (0) (r^|r 0 ,/ 0 ) 
47i 4x jUp 

< 

G, (1) (r^|r 0 ^ 0 ) = -^G (1) (r,/|r 0 ^ 0 )--^-— s-VG w (r,t\r 0 ,t 0 ) 
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where \x tr is the transport scattering property defined as 

jx fr = fi a + f/ s9 and /J s = (1 - g) // 5 with g the mean cosine of the forward scattering 

angle. It can be shown from Eq. (11) that the first order correction to the Green's 
function is: 

5 (17) 
G (1) (r,^ 

+ -VvG (0) (r 3 /|r'/).^G (0) (r' 9 ^|r 0? g 

3 Mtr 



The above equation is identical to the correction calculated directly from the diffusion 
equation with perturbation theory. 

The present invention involves the series expansion method and the algebraic 

10 reconstruction technique (ART) for optical CT. The volume to be imaged is split into 
N x N x N cubic voxels. Each voxel is assigned a value representing the local average 
of the absorption distribution. The imaging problem is 2D in the case of X-ray. The 
linear attenuation in Eq. (2) is then simplified to a summation of the voxel values along 
the source-detector line. The ray sum can be exactly expressed as: 

15 (18) 

N 2 

y J = i:b ij w iJ ii r 

In Eq. (18), y } is the data of the yth measurement: w is a geometric factor related 
to the oblique angle of the y*th measurement, and it equals the segment length of the y'th 
ray within voxel z; and \i { is the local average of the absorption within voxel L The 
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10 



summation on the right hand side of Eq. (18) is that of the voxels on the detection 
plane. In the case of X-rays, the factor b y is given as: 



As illustrated in Figure 2A, through the introduction of the weighting factor b lJ9 the plane 
summation in Eq. (18) is actually a line summation. Those voxels which fall on the line 
give 100% contribution to the X-ray function, while those which fall away from the line 
give 0% contribution. 

Our extension to optical CT is through the weighting factor b iJt . Because of 
scattering photons have different probabilities to follow different paths. Therefore, the 
contribution of each voxel to a source-detector measurement is weighted by the density 
of photon paths across it. The value of b for the optical case is the value of the photon 
path probability (Figure 2B): Eq. (18) can then be directly applied. Recalling that the 
photon path is a 3D tube, summation over three dimensions is required. The generalized 
equation is: 



(19) 




1, the j th ray crosses pixel z; 
0, other. 



(20) 




This is exactly the discrete form of Eq. (14) at a fixed time i, if we set 



(21) 
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' yj = -[G f (r y ,/|r 0 ^ 0 )-G{ 0) (r y ,f|r 0 ^)], 

^ PyWy = PSF( rj , t; r, ; r 0 r,0) • Av, 

with Av the volume element of each voxel. Equation (20) can be rewritten in a compact 
matrix form 

(22) 

5 y = Rx + n, 

where y is the measurement vector (dimension M), x is the image vector (dimension N 3 ), 
R is the projection matrix containing geometric and photon path information (dimension 
M x N 3 ), and n denotes the error vector. 

In image reconstruction, apply the additive ART of X-ray CT directly to optical 

10 CT. Since the imaging problem can be either overdetermined or underdetermined, it is 
inappropriate to solve the equation y = Rx directly, as it may not have any solution at all, 
or it may have many solutions but none is appropriate. The estimation of the image 
vector is usually performed using optimization criteria, based on the error between 
forward model predictions and experimental data, as well as a priori information about 

15 the imaging region. One example is the so-called regularized least-squares criterion, 
which is a special case of the Bayesian estimate. The reconstruction of the sought-after 
image is performed through minimizing the function: 

(23) 

r 2 ||y-Rx|| 2 +||x-5jx 0 | 2 . 

20 In Eq. (23), r is a parameter often called the signal-to-noise ratio in the literature 

8|i 0 is the pre-knowledge for the image vector (N 3 x 1) and also used as the initial 
estimate of the image, and il.ji 2 denotes the module square for the vector. Keeping the 
second term small ensures that the picture is not too far from the pre-knowledge, and 
keeping the first term small ensures that the picture is consistent with the measurements. 
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One iterative procedure to minimize Eq. (20) is to introduce two sequences of 
vectors, x (k) and u (k) , of dimensions A? and M 9 respectively. Initially, x (0) is set equal to 
6\x 0 and u (0) to a zero vector. The iterative step (G.T. Herman, Image Reconstruction 
From Projections: The Fundamentals of Computerized Tomography (Academic Press, 
5 New York, New York, 1980)) which is incorporated herein by reference in its entirety is 
given by: 

(24) 

r x (w) =x (*) +rc w R 

j k 

u (* + i) =u w +c (*) e . ' 



where 



10 



(25) 



„(*) 



= A.- 



-u, 



1 + r 2 R. 



with i k = k{mod JV 3 ) + 1 . if, the transpose of the z-th row of the matrix R. e, the M 
dimensional column vector with 1 in the z-th row and 0 elsewhere, and X is a real 
number called relaxation parameter. It has been proven that as long as 0<l<2 the 

15 sequence {x (k) } converges in the limit to the minimizer of Eq. (23). 

In order to demonstrate this optical CT technique, measurements were made in a 
cubic glass container 6.35 cm on a side, filled with a tissue-like scattering medium. The 
cubic geometry was selected to simplify the theoretical modeling because of its regular 
boundaries. The scattering medium was a stock solution prepared with 1.072 |im 

20 diameter polystyrene beads (PolyScience, Inc.) at 0.27% concentration. The scattering 
and adsorption properties of the medium were determined in real time by fitting the 
transmitted diffuse light. A schematic diagram of the system 100 is presented in Figures 
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3 A and 3B. The light source 102 was a Coherent Mira 900 mode-locked Ti:sapphire 
laser operated in femto-mode and pumped by a Coherent Innova 400 multiline argon ion 
laser 104. The wavelength was 800 nm, and the repetition rate was 76 MHZ. The pulse 
width was -150 fs. The detection system 106 was a Hamamatsu streak camera C5680 
5 with M5675 Synchroscan Unit. A small portion of the laser beam was deflected by a 
quartz plate to a fast photodiode 108 (Hamamatsu CI 808-02), which generates the 
triggering signal 1 10 for the streak camera. The cubic glass container was mounted on a 
translation stage 112 used to actuate relative movement between the source, the detector 
and the object to be scanned. A Pentium-H computer 120 served as a central controlling 

1 0 and data acquisition unit. In order to monitor the laser power drift during the scan, a 
DT2801-A card (Data Translation, Inc.) was installed on the computer and programmed 
to record the laser power voltage from the control box of the light source, and it also 
served to drive the stepper motor of the translation stage. Total automation of data 
acquisition for a full line scan was achieved through programming the user interface of 

15 the streak camera software. The streak camera was operated in analog mode. It 
converted the temporal evolution of light signals into vertical streak images. The 
transmitted light was collected with four coherent fiber bundles 130. Each fiber bundle 
had a 500 |im core diameter and consisted of ten thousand single mode silica fibers. The 
proximal ends of the four fibers were bundled together to increase the collection area. 

20 The overall time resolution of the detection system was set to 30 ps. 

The container was placed in the sample holder on the translation stage. Multiple 
objects were embedded inside the turbid medium at fixed positions. The sample holder 
was designed so that the top, bottom, left and right boundaries of the container were 
totally black. The laser pulses were delivered into the medium at the front side, and the 

25 transmitted signals were collected at the opposite side. The incoming laser beam and the 
collection fibers were aligned in a coaxial geometry. Surface scans were conducted on 
the XZ and YZ directions, so that a 3D image of the absorption distribution could be 
reconstructed. The scan area was a 4cm x 4cm square along each direction. The scans 
were 2 mm per step in the horizontal and vertical directions, resulting in 2 projections, 
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882 total scan measurements. The data acquisition time for each point was 8s. Either 
one or two opaque objects were embedded in the medium. In each case, the scans were 
first carried out for the XZ plane, and then the container was rotated by 90° and the scans 
continued for the YZ plane. 
5 The calculation of the point spread function requires knowledge about the 

average scattering and absorption properties of the scattering medium and ju a ) . 

To determine \x ' s and ju a , the coaxial transmission signal of a 2.2 ns time window was 

collected, with the absorbers removed from the scattering medium and the source- 
detector aligned at the center of the X-surface. This time-dependent curve was fitted 
10 using the diffusion approximation solution with zero-boundary condition. The best fit 

was given by jj,^ = 738cm' 1 and fi a = 0.03cm" 1 . These values are similar to those of 

human breast tissue. 

In these measurements, data was collected for two configurations. In the first 
case, an 8 mm diameter black sphere was placed at the center of the container ("one 

1 5 object configuration"). In the second case, a black sphere 8 mm in diameter was placed 
at (0.56, 56.-0.56) cm from the center, and a second black sphere 6 mm in diameter was 
placed at (-0.56.-0.56.56) cm from the center ("two object configuration"). For each 
case, scanning was performed on the surfaces in 2 mm steps, and 882 time evolution 
curves of the transmission signals were measured, one for each scanning point. Because 

20 the source-detector distance remained the same for all the 882 measurements, the time 
scale was the same and the intensities can be compared at the same time point. When 
the source-detector line crosses the absorber there is a drop in signal level, and vice 
versa. The projection intensity diagrams can be easily created by plotting the contour 
map of the intensities of all the 882 time evolution curves at the same time point. To 

25 achieve higher counts and reduce noise, such intensities can be summations of counts 
over a time window comparable to the time resolution of the detecting system. The 
selection of the time point is based on the trade-off of spatial resolution and 
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the absorbers have weaker intensity. The absorbers appear as shadows on the projection 
intensity diagrams. The projections of one embedded absorber (Figures 5A and 5B) are 
different from those of two embedded absorbers (Figures 5C and 5D). Figures 5A-5D 
are referred to as the zero-th order images. They already show features such as the 
5 central positions of the embedded objects, though the images are scrambled due to 
scattering. 

Before processing the data of Figures 5A-5D, artifacts at the edges due to tiny air 
bubbles and asymmetry of the surface were removed manually, and a grid was created. 
In this analysis choose the 2-mm scanning step as the size of each grid, resulting in 

10 21 3 =9261 voxels. The inverse is underdetermined, as 9261 voxel values are 
reconstructed from 882 measurements. 

The image reconstruction procedure was then applied in two steps. First, the 
straight line PSF of X-ray CT was applied to the data using ART (Eqs. (24) and (25)) 
with an initial estimate of 5|i = 0 (See Eq. (23)) and "signal-to-noise" r = 60. Second, 

15 the resulting image, obtained with the straight line PSF, was then used as an initial 
estimate and applied to the data using a particular photon migration PSF. The 
"signal-to-noise ratio" in Eq. (23) was again set to r = 60, and the relaxation parameter in 
Eq. (25) was set to X = 1. The reconstructions were computed with point spread 
functions calculated from both the diffusion approximation and the causality corrected 

20 solutions (see below). The number of iteration steps was set to 2 x 10 5 . On a Pentium-II 
450MHz PC, each reconstruction took -40 seconds. The number of iteration steps was 
increased to 2 x 10 6 without observing much improvement. 

The second step involves the calculation of the PSF from solutions to the 
transport equation. The diffusion approximation solution is widely used in the literature. 

25 However, it is well known that this solution violates causality and breaks down in the 
early time regime. Solutions incorporating causality have been worked out for models 
based on random walk and path integral theories. A Green's function which 
incorporates causality and is valid for early arriving photons has been used. This 
Green's function is constructed from the diffusion approximation Green's function 
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The second step involves the calculation of the PSF from solutions to the 
transport equation. The diffusion approximation solution is widely used in the literature. 
However, it is well known that this solution violates causality and breaks down in the 
early time regime. Solutions incorporating causality have been worked out for models 
5 based on random walk and path integral theories. A Green's function which 
incorporates causality and is valid for early arriving photons has been used. This 
Green's function is constructed from the diffusion approximation Green's function 

G(r, *|r 0 , t 0 ) by replacing the time t 0 at which the pulse is injected into the medium with 

t 0 +t iT with t iT the time of flight for unscattered photons to propagate from the source to the 

10 detector. Physically, this procedure takes into account the fact that diffusion starts only 
after the light pulse traverses the medium. 

Figures 6A-6D show the point spread functions calculated with the diffusion 
approximation solution and with the causality corrected solution. Note that all contour 
maps in Figures 6A-6D correspond to the time window of the experimental data for the 

15 two-object case. The point spread function calculated using the diffusion approximation 
is wider than that using the causality modified procedure. 

The image reconstruction results for X-ray (straight line PSF), diffusion 
approximation, causality correction, and the actual configuration are presented in the 
contour plots of Figures 7A-7D and Figures 8A-8D. Figures 9A-9D present 

20 characteristic slices of Figure 8C in 2D contour maps. The slices are picked at planes z 
= -10 mm z = -6 mm z = 2mm and z = 10 mm. As clearly seen in Figures 7A-7D and 
8A-8D, the images reconstructed with the straight line path (X-ray CT) do not have high 
resolution, while the images reconstructed with the diffusion approximation are too 
small compared with the actual size of the spheres, especially for the two object case, for 

25 which the smaller object is invisible. Due to the noise in the experimental data, the 
reconstructed images for the one object configuration exhibit some distortions. 
Remarkably, for the two object configuration the image reconstructed with the causality 
correction give correct sizes of the embedded objects. As shown in Figure 9, the voxel 
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values off the objects are nearly zero which naturally results from the inverse procedure. 
In Figure 8C, the size and the position of the 8-mm object and the size of the 6-mm 
object are correct, while the position of the 6-mm object is off from its actual position by 
2 mm. This may indicate that a cross-talk exists in the inverse when two objects are 
5 embedded. 

The point spread function calculated with the conventional diffusion 
approximation solution is too wide for the early time window of our data. The 
reconstructed images in this case are too small compared with the actual size of the 
objects. Especially for the two object configuration, the smaller absorber is basically 

10 invisible in the reconstructed image. On the other hand, the same calculations done with 
the causality correction provide improved resolution. 

Figure 10 illustrates a process sequence 200 in accordance with a preferred 
embodiment of the invention. After storing 202 a light distribution function and/or any 
reference data in electronic memory and programming 204 the computer to process the 

1 5 collected image data for a particular type of anatomy or tissue structure such as a 
concerous lesion, the patient or biopsied sample is scanned 206 with an endoscope or 
probe. The collected data is processed 208 to generate an image for display and for 
further processing 210. 

While this invention has been particularly shown and described with references 

20 to preferred embodiments thereof, it will be understood by those skilled in the art that 
various changes in form and details maybe made therein without departing from the 
scope of the invention encompassed by the appended claims. 



